#!/usr/bin/env python
# coding: utf-8

"""

作者：张陈亚
微信公众号名称：张陈亚
微信号: zy10178083

【淘宝店铺地址】：https://shop63939877.taobao.com/
【淘宝店铺名称】：张陈亚

【项目实战合集导航】：
https://docs.qq.com/sheet/DTVd0Y2NNQUlWcmd6?tab=BB08J2


关于项目代码运行相关问题的答疑.txt 这个说明一定要看一下，避免运行报错。

"""

# 导入需要的包
import pandas as pd
import numpy as np
from arch import arch_model
from scipy.stats import shapiro
from scipy.stats import probplot
from statsmodels.stats.diagnostic import het_arch
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.stats.diagnostic import acorr_ljungbox
from matplotlib import pyplot as plt
from statsmodels.stats.diagnostic import acorr_ljungbox

plt.style.use('fivethirtyeight')
plt.rcParams['xtick.labelsize'] = 10
plt.rcParams['ytick.labelsize'] = 10

# 读取数据
df = pd.read_csv('F:/LTFX/600059.CSV', index_col='date')
# 用Pandas工具查看数据
print(df.head())
# 查看数据集摘要
print(df.info())

df = df.drop(columns=['open', 'high', 'low', 'volume'])

df['pct_change'] = 100 * df['close'].pct_change()
df.dropna(inplace=True)
# 股票收盘价趋势分析
df['close'].plot(figsize=(10, 5), title=f'TDG Closing Price 2013-2018')
plt.show()
# 股票每日收益率趋势分析
df['pct_change'].plot(figsize=(10, 5), title=f'TDG Percent Change in Closing Price')
plt.show()
# 股票每日收益率自相关图
acf = plot_acf(df['pct_change'], lags=30, title=f'TDG Percent Change Autocorrelation')
plt.show()
# 股票每日收益率偏自相关图
pacf = plot_pacf(df['pct_change'], lags=30, title=f'TDG Percent Change Partial Autocorrelation')
plt.show()

# 股票每日收益率白噪声检验
#ljung_res = acorr_ljungbox(df['pct_change'], lags=40, boxpierce=True)
#print(f'Ljung-Box p-values: {ljung_res[1]}')
#print(f'Box-Pierce p-values: {ljung_res[3]}')
lb_stat, lb_pvalue, bp_stat, bp_pvalue = acorr_ljungbox(
    df['pct_change'], lags=40, boxpierce=True)

# 打印结果
print(f'Ljung-Box p-values: {lb_pvalue}')
print(f'Box-Pierce p-values: {bp_pvalue}')

def ts_plot(residuals, stan_residuals, lags=50):
    residuals.plot(title='GARCH Residuals', figsize=(15, 10))
    plt.show()
    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(20, 10))
    ax[0].set_title('GARCH Standardized Residuals KDE')
    ax[1].set_title('GARCH Standardized Resduals Probability Plot')
    residuals.plot(kind='kde', ax=ax[0])
    probplot(stan_residuals, dist='norm', plot=ax[1])
    plt.show()
    plot_acf(stan_residuals, lags=lags, title='GARCH Model Standardized Residual Autocorrelation')
    plt.show()
    plot_pacf(stan_residuals, lags=lags, title='GARCH Model Standardized Residual Partial Autocorrelation')
    plt.show()


# 建立
garch = arch_model(df['pct_change'], vol='ARCH', p=1, q=1, dist='normal')  # vol='GARCH' 改成 vol='ARCH'
fgarch = garch.fit(disp='off')
resid = fgarch.resid
st_resid = np.divide(resid, fgarch.conditional_volatility)
ts_plot(resid, st_resid)
print('****************fgarch.summary()*********************')
print(fgarch.summary())

arch_test = het_arch(resid, nlags=50)
shapiro_test = shapiro(st_resid)

print(f'Lagrange mulitplier p-value: {arch_test[1]}')
print(f'F test p-value: {arch_test[3]}')
print(f'Shapiro-Wilks p-value: {shapiro_test[1]}')


# 模型优化
def gridsearch(data, p_rng, q_rng):
    top_score, top_results = float('inf'), None
    top_models = []
    for p in range(len(p_rng)):
        for q in range(len(q_rng)):
            model = arch_model(data, vol='GARCH', p=p_rng[p], q=q_rng[q], dist='normal')
            model_fit = model.fit(disp='off')
            resid = model_fit.resid
            st_resid = np.divide(resid, model_fit.conditional_volatility)
            results = evaluate_model(resid, st_resid)
            results['AIC'] = model_fit.aic
            results['params']['p'] = p_rng[p]
            results['params']['q'] = q_rng[q]
            if results['AIC'] < top_score:
                top_score = results['AIC']
                top_results = results

            elif results['LM_pvalue'][1] is False:
                top_models.append(results)

    top_models.append(top_results)
    return top_models


# 模型评估
def evaluate_model(residuals, st_residuals, lags=50):
    results = {
        'LM_pvalue': None,
        'F_pvalue': None,
        'SW_pvalue': None,
        'AIC': None,
        'params': {'p': None, 'q': None}
    }

    arch_test = het_arch(residuals, nlags=lags)
    shap_test = shapiro(st_residuals)

    results['LM_pvalue'] = [arch_test[1], arch_test[1] < 0.05]
    results['F_pvalue'] = [arch_test[3], arch_test[3] < 0.05]
    results['SW_pvalue'] = [shap_test[1], shap_test[1] < 0.05]

    return results


p_rng = list(range(1, 3))  # 为了缩短运行时间 我把31 改为3
q_rng = list(range(1, 5)) # 为了缩短运行时间 我把41 改为5
df['dif_pct_change'] = df['pct_change'].diff()  # 一阶差分 序列变为平稳非白噪声序列
df = df.reset_index()
top_models = gridsearch(df['dif_pct_change'].iloc[1:, ], p_rng, q_rng)
print('*****************top_models*******************')
print(top_models)

# 建模：使用最优的 p  q 值
garch = arch_model(df['pct_change'], vol='GARCH', p=17, q=25, dist='normal')
fgarch = garch.fit(disp='off')
resid = fgarch.resid
st_resid = np.divide(resid, fgarch.conditional_volatility)
ts_plot(resid, st_resid)
arch_test = het_arch(resid, nlags=50)
shapiro_test = shapiro(st_resid)
print(f'Lagrange mulitplier p-value: {arch_test[1]}')
print(f'F test p-value: {arch_test[3]}')
print(f'Shapiro-Wilks p-value: {shapiro_test[1]}')
print('*****************fgarch.summary()*******************')
print(fgarch.summary())

print('预测未来5天，其预测方差：\n')
forecasts=fgarch.forecast(horizon=5)
print(forecasts.variance.dropna().head())